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Abstract 

We are interested in nonlinear hyperbolic systems in nonconservative form arising in fluid dynam- 
ics, and, for solutions containing shock waves, we investigate the convergence of finite difference 
schemes applied to such systems. According to Dal Maso, LeFloch, and Murat's theory, a shock 
wave theory for a given nonconservative system requires prescribing a priori a family of paths 
in the phase space. In the present paper, we consider schemes that are formally consistent with 
a given family of paths, and we investigate their limiting behavior as the mesh is refined, we 
first generalize to systems a property established earlier by Hou and LeFloch for scalar conser- 
vation laws, and we prove that nonconservative schemes generate, at the level of the limiting 
hyperbolic system, an convergence error source-term which, provided the total variation of the 
approximations remains uniformly bounded, is a locally bounded measure. This convergence 
error measure is supported on the shock trajectories and, as we demonstrate here, is usually 
"small" . In the special case that the scheme converges in the sense of graphs — a rather strong 
convergence property often violated in practice — then this measure source-term vanishes. We 
also discuss the role of the equivalent equation associated with a difference scheme; here, the 
distinction between scalar equations and systems appears most clearly since, for systems, the 
equivalent equation of a scheme that is formally path-consistent depends upon the prescribed 
family of paths. The core of this paper is devoted to investigate numerically the approximation 
of several (simplified or full) hyperbolic models arising in fluid dynamics. This leads us to the 
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conclusion that for systems having nonconservative products associated with linearly degenerate 
characteristic fields, the convergence error vanishes. For more general models, this measure is 
evaluated very accurately, especially by plotting the shock curves associated with each scheme 
under consideration; as we demonstrate, plotting the shock curves provide a convenient approach 
for evaluating the range of validity of a given scheme. 

Key words: nonconservative hyperbolic system, shock wave, family of paths, equivalent equation, 
convergence error measure, formally path-consistent scheme. 
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1. Introduction 

A number of non-conservative hyperbolic models have been introduced in fluid dynam- 
ics to serve as (simplified) models of two-phase or two-layer flows. Our objective in the 
present paper is to address the fundamental question whether finite difference schemes 
for nonconservative systems converge toward correct weak solutions containing shock 
waves. Addressing this important issue requires detailed numerical computations which 
we carry out here. The nonconservative hyperbolic systems under consideration have the 
general form 

u t + A(u) u x = 0, u — u(t, x) e R N , (1.1) 

where it is the vector-unknown and A — A(u) is a smooth, N x N matrix-valued map A 
which admits real eigenvalues \\ < . . . < Ajy and a basis of eigenvectors r±, . . . , rjy. We 
are interested in solving the initial value problem associated with some initial condition 

u(0,x) = u (x), i£R. (1.2) 

The solutions of nonlinear hyperbolic systems arc generally discontinuous; due to the 
non-divergence form of the equations the notion of solutions in the sense of distributions 
can not be used, and weak solutions satisfying (1.1) are defined in the sense introduced 
by LeFloch [18,19,20,21,23] and Dal Maso, LeFloch, and Murat [13] (cf. Section 2 below 
for a brief review of the theory) . 

Generally speaking, solutions to (1.1) depend upon regularization mechanisms; for 
instance, different approximation schemes may converge toward different solutions, and 
for this reason in developing the well-posedness theory, higher-order regularization effects 
such as viscosity, capillarity, relaxation terms, must be taken into account in the modeling. 
For instance, for continuous models one may consider the regularization 
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ul + A(u^ul = R e , 



(1.3) 



where R e depends upon higher-derivatives of u together with (one or several) small-scale 
parameter(s) e; the physically meaningful solutions are defined as the singular limits 



Furthermore, as established in [19], shock waves in such solutions are determined by 
traveling wave solutions to (1.3), that is, Rankine-Hugoniot relations for shock waves are 
determined from the given regularization. 

In the present paper, we demonstrate that while, for certain simplified models, solu- 
tions are actually stable upon regularization and the detailed knowledge of the right-hand 
side of (1.3) is unnecessary, however for general systems such as the "full" systems of 
two-phase flows, the general DLM theory is necessary. Still, as pointed out by Hou and 
LcFloch [16] — who focused attention on the same issues for nonconservative formula- 
tions of scalar hyperbolic equations — and by Hayes and LeFloch [14,15] and LeFloch 
and Mohammadian [25] — who studied the effect of diffusive and dispersive terms — the 
effects of the regularization R e may be difficult to pinpoint in practice. In view of the 
fact that the models under study are derived from modeling approximation assumptions, 
this fully justifies the use of a numerical strategy based on a direct discretization of the 
nonconservative hyperbolic models (1.1). Our conclusions, therefore, justify to search for 
robust and efficient high-order schemes for the approximation of nonconservative systems. 
In particular, Berthon and Coquel [1,2] and Chalons and Coquel [11], have introduced 
various numerical strategies for models of complex fluid flows including turbulence mod- 
els, while Pares [30] and Munoz-Ruiz and Pares [28] have developed many important 
applications. 

Finally, for anther standpoint to the theory and the numerical analysis of nonconser- 
vative products, we refer to Berthon, Coquel, and LeFloch [3] who connected the theory 
of nonconservative products with the concept of a kinetic relation [22] . They introduced 
a general framework to handle nonconservative systems; this framework encompasses 
a large number of examples arising in the applications. In particular, they rigorously 
analyzed a typical model of turbulent fluid dynamics by establishing the existence and 
properties of a physically relevant family of traveling waves and deriving the correspond- 
ing kinetic function. 

2. The convergence error measure 

2.1. DLM familes of paths and nonconservative products 

Let f2 be an open subset of and g : f2 — > be a smooth mapping. Given a function 
with bounded variation u : R — > f2, the definition introduced by Dal Maso, LeFloch, 
and Murat [13] allows one to define products of the form g(u) ^ provided a family of 
Lipschitz continuous paths $ : [0, 1] x f2 x f2 is prescribed, which must satisfy certain 
natural regularity conditions, in particular 



u := lim u e . 

£^0 



$(0;ui,u r ) = ui, 



${l;ui, u r ) = u r , 



(2.1) 



$(s;ui,ui) = m 



(2.2) 
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for all s € [0, 1] and m, u r e ft. The nonconservative product, denoted by [g(u) ^] $ , is 
defined as a bounded measure, which is absolutely continuous with respect to the total 
variation measure of the function u and, in particular, coincides with the distributional 
derivative -£-f in the special case of a conservative product, for which g(u) = Df(u) for 
some /. 

The DLM theory was applied to nonconservative systems of the form (1.1); the Rie- 
mann problem was solved and, later, the general Cauchy problem [24]. In the course of 
this analysis, the notion of ^-completion (X, U$) of the graph of a BV function u was 
introduced. The key stability result in [13] is the following: if u A : K — > ft is a sequence 
of BV functions with uniformly bounded amplitude and total variation 

S up\u A \ + TV R (u A )<l, (2.3) 

R 

converging almost everywhere to a limit function u : M — > ft, then a sufficient condition 
for the corresponding nonconservative products to converge 

ff ( U A )^l -[«?(«) ^1 (2.4) 

in the weak-star sense of measures, is that their ^-completions (^ A , t/^) converge in the 
uniform distance of graphs precisely to the ^-completion [X, L/$) of the limit u. 

2.2. ^4 c/ass of finite difference schemes 

For the approximation of nonconservative systems we introduce here a general family 
of numerical schemes which includes, in particular, three classes of schemes of particular 
interest: Godunov, Roe, and Lax-Friedrichs. 

For the discretization of the initial value problem (1 . 1 )- (1 .2) , we introduce computing 
cells Ii = [xj_i/2,Xj+i/2] and, for simplicity in the presentation, we assume that these 
cells have constant size A = Ax. We also define x i+ i/ 2 = iAx and Xj = (i — 1/2) Ax, 
the latter being the center of the cell l im Finally we denote by At the (constant) time 
length and we set t n = nAt. We denote by u" the approximation of the cell averages of 
the exact solution provided by the numerical scheme: 

< - — / u(t n ,x)dx. 
We are interested in schemes of the general form 

< +1 = <-^ x {K^ + Ku,^ (2-5) 

where 

ACv2 = Af± («?- 9 .-".«?+ P )- 
From now on a DLM family of paths $ for the nonconservative system (1.1) is fixed. 
Following Pares [30], we consider formally path-consistent schemes, i.e. schemes that are 
consistent with the family of paths <& in the following sense: M~ and M + are Lipschitz 
continuous mappings from £l p+q+1 to il satisfying: 

M ± (u,...,u) = 0, it eft, (2.6) 



5 



and for every m £ CI, i — —q, . . . ,p, 

M ■ • ■ , u p ) + M + (u_ q , . . . , u p ) = J A($(s;jt ,u 1 )) — (s;u ,Ui)ds. (2.7) 

These conditions provide a generalization to the concept of conservative scheme intro- 
duced by Lax for systems of conservation laws and, for this reason, were originally called 
"path-conservative" in [30]. As we will see, this definition -although quite natural- needs 
to be handled carefully. 

It is convenient to assume some particular structure on the given family of paths. 
Precisely, we assume that the matrix A and the paths <& satisfy the following restrictions: 
(Rl) Given an integral curve 7 of a linearly degenerate field and ui,u r £ 7, the path 
<i>(s; m, u r ) is a parametrization of the arc of 7 connecting m and u r . 

(R2) Given an integral curve 7 of a genuinely nonlinear field and ui,u r £ 7, with X(ui) < 
X(u r ), being X(u) the corresponding eigenvalue, the path <I>(s;u ; ,u r ) is a parame- 
trization of the arc of 7 connecting ui and u r . 

(R3) Let us denote by WJ 3 cClxfl the set of pairs (ui,u r ) for which the Riemann problem 

u t + A(u)u x = 0, 

ui, x < 0, ( 2 - 8 ) 



u(x, 0) 



x > 0, 



has a unique self-similar weak solution composed by N (possibly trivial) simple waves 
connecting N + 1 intermediate constant states 

u = ui, m,... ,Un-i, u n = u r . 

Given (ui,u r ) £ 3W, the curve described by the path <&(•; m, u r ) is equal to the union 
of those corresponding to the paths <&(•; Mj-i, Uj), j — 1, . . . , N. 

We now present several schemes of particular interest in the present paper. First of all, 
the Godunov method for the nonconservative system (1.1) takes the form (2.5) with 

KiT/2 = [ A ^ «?. <+i/ 2 )) ^ «?. <+i/ 2 ) da, 

K+1/2 =J Q ^(*(*»«"+i/2.«?+i))^-(*;«"+i/2.«"+i)d». 

where u™ +1 ^ 2 is the (constant) value at x = of the solution of the Riemann problem 
consisting of (1.1) with initial condition: 



u(x, 0) = 



*i 1 



x < 0, 



x>0. 



In the derivation it is important to assume a "CFL-1/2 condition", as noted in [28]. 

In the case in which the solution of the Riemann problem is discontinuous at x = 0, 
such discontinuity has to be stationary and we can replace u™ +1 / 2 cither by the limit of 
the solution to the left or to the right of 0. 
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Secondly, Roe methods provide linear approximate Riemann solvers. In view of the 
framework in [13] the following generalization to Roe's standard approach was proposed 
by LeFloch [21]. Given a family of paths <f>, a function A$ : il x il i— > Mat x at(]R) is called 
a (generalized) Roc linearization if for all ui,u r e ft the following properties hold: 

(i) A$(m,u r ) has N distinct real eigenvalues, 

(ii) A<s>(ui,ui) = A(u t ), 

f 1 <9$ 

(iii) A<s>(ui,u r ) ■ (u r - ui) = / A($(s; u u u r )) — (s; u u u r ) ds. 

Jo ds 

Once a Roe linearization has been chosen, the corresponding Roe scheme takes the 
form (2.5) with 



M 



+ 1/2 ~~ /1 j+l/2 \ a i+l 



A/f™' + — A n ' + ■ (ii n — n n \ 
1V1 i+l/2 ~ A i+l/2 \ U i+l 



where A? +1/2 = Aj,«,< +1 ), 



and 



n,± 
i+1/2 



A i+l/2 



(K+l/2,l)~ 



/2,N 



) ± 



• n,± 



+ 1/2^+1/2 



(^+1/2) 



Here, £<™ +1 / 2 denotes the diagonal matrix whose coefficients are the eigenvalues of A™ +1 ^ 2 , 

J+l/2,1 ^ A i+l/2 : 2 ^ 

associated eigenvectors. 



< A™ +1 ^ 2 N , and 3C" + i/2 & N x N matrix whose columns are 



Third, a generalization of the classical Lax-Fricdrichs method to (1.1) is given by (2.5) 
with the choice: 



Mf+} /2 - jT' K, < +1 )) ^ (*; «?, <+i) ds, 



(2.9) 



where 



A-( U ) = 



Ax 
'At 



Id + A(u) 



A + (u) 



Ax 
~At 



Id + A{u) 



being Id the N x N identity matrix. 

2.3. Convergence to a nonconservative system with measure-source term 

We denote by u A the sequence of piecewise constant approximate solutions generated 
by a finite difference scheme of the form described above. By extending the arguments 
in Hou and LeFloch [16] we can prove: 

Claim 1 Consider a nonconservative hyperbolic system (1.1) together with a given fam- 
ily of paths <&. Suppose that u A is a sequence of approximate solutions constructed by 
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one of the finite difference schemes described in Subsection 2.2 and satisfying the bounds 
(2.3) uniformly in time. Suppose that the scheme is formally consistent with the family 
of paths Then, given any subsequence of u A converging almost everywhere to some 
limit, denoted by v, the following holds: 
(i) There exists a bounded measure \i v : R + x R — > R N (called the convergence error 
measure,) such that the limit v satisfies the following hyperbolic system with source- 
term 

Vt+[A(v)v x ] 9 = ii v . (2.10) 

(ii) Moreover, when the ^-completion of the graphs of u A converges in the uniform 
sense of graphs towards the ^-completion of the limit v, i.e. 

(X A ,Ui)^(Y,V*), 

then the convergence error measure [i v vanishes identically and v is a weak solution 
to the system 

vt+[A(v)v x ] 9 = 0. (2.11) 



The above result can be interpreted as a "nonconservative extension" to the classical 
Lax-Wendroff theorem for systems of conservation laws [17]. It should be observed that 
the convergence in the sense of graphs is very strong; it does hold for the Glimm and 
front tracking schemes, but usually fails for finite difference schemes. 

Our main objective in the present paper is to investigate the source and the amplitude 
of this convergence error, which can be measured in terms of the measure fi v or, equiva- 
lcntly, in terms of the Rankine-Hugoniot curves associated with the given scheme. Indeed, 
computing numerically the shock curves associated with various schemes of interest is 
one of the main purposes of this work. 

Proof. We follow Hou and LeFloch [16] and decompose the scheme into a part that 
converges to the hyperbolic system and an error term, and we then rely on stability 
results established in the general DLM theory [13]. We need to prove that 



oo poo 



v(t,x)(p t (t,x)dxdt- ([A(v(t,-))v x (t,-)] 9 ,<p(t,-))dt = (2.12) 

J-oo JO 

for all compactly supported test-function p = <p(t, x). 
We set 

<p? = <p(t n ,Xi), ip? +1/2 = <p{t n ,x i+1/2 ) 

and we multiply the discrete equation (2.5) by y>™. After summing over i and n, and then 
applying summation by parts we obtain the identity 



oo oo n „_1 



At 

(2.13) 



OO OO 



n— i- 



We want to prove that equation (2.12) can be obtained by passing to the limit as A 
tends to in (2.13). The convergence of the first term in (2.13) to the corresponding one 
in (2.12) is obtained as in the conservative case. Concerning the second term, since M ± 
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are Lipschitz continuous and the states w™ are uniformly bounded, we only have to study 
the convergence of 



oo oo 



n— i— — oo 

oo oo , „i v 
n=0i=-oo ' 

where we have used our assumption that the scheme is formally consistent with the given 
family of paths <&. 

Using the definition of the product A(u A (t, -))u A (t, •) as a measure this term can be 
rewritten as follows: 

/•OO 

/ <[4(u A (V))«£(t, ■)]*,¥>(*,■)><** 
Jo 

+ E E / / A($( a; <,«? +1 ))$ 8 (*;<,«? +1 )d* )(^ 1/2 -<p(t,x i+1/2 ))dt. 

The second summand trivially converges to 0. The proof is concluded by using the main 
stability result and error estimate in [13]. Indeed, if the sequence u A converges in the 
sense of graphs, then, as recalled in (2.4) the nonconservative product converges in the 
weak-star sense of measures, 

/ ([A(u A (t,-))u A (t,-)]^v(t,-))dt^ / ([A(v(t,-))v x (t,-)]*,v(t,-))dt 
Jo Jo 

and we recover the exact system (2.11). In the general case, the <I>-graph completion 

(X A , U A ) of u A converges to some limiting graph (X, U), 

(X A ,U A )^(X,U), 

which is such that its projected BV function coincides with the pointwise limit u of u A . 
In turn the limit nonconservative product is based on the BV function v and we obtain 



Jo 



([A(u A (t,-))u A (t,-)]^^-))dt 

A(U(t, -))U(t, •)) dt= I ( [A(v(t, -))v x (t, •)] 4 + » v (t),<p(t, •)) dt 



for some time-dependent measure /j, v . This completes the proof. □ 

3. Equivalent equations for nonconservative systems 

3.1. Derivation of the equivalent equation 

In this section, we derive the equivalent equations corresponding to the Lax-Friedrichs- 
typc scheme introduced in Subsection 2.2, which reads 

^ (V 1 - 2 K-i + <+0) + aL(jf ^ (*(*;«?_!,«?)) ^(s;u?-i,K)ds 

+ J\ ($(*; <, < +1 )) g( S ; <, < +1 ) ds) = 0. 



We consider first the equivalent equation at second order. Performing a formal Taylor 
expansion in the scheme we obtain: 

At ( Ax 2 Ax \ 

v t + A(v)v x + — lv tt - v xx + — I 2 {v) J = 0, 

h{v)= [ DA(v)(D ul $-v x ,D ul $ s -v x ) ds+ [ DA (v) (D Ur & ■ v x , D Ur $ s ■ v x ) ds, 
Jo Jo 

where the following notation has been used: given two vectors v = [v i, . . . , vn] t , w = 
[wi, . . . , wn] t , DA(u)(v, w) represents the derivative of A(u) in the direction of the vector 
v (which is a matrix) applied to the vector w, i.e. DA(u)(v, w) — (^2^=1 v kd Uk A(u)j ■ w, 
where d Uk A(u) is the N x N matrix whose entry is d Uk a,ij(u). 

As usual, we write the ^derivatives using x-derivatives in order to obtain a new mod- 
ified equation. Using the equality 

(A 2 (v)v x ) x = DA(v) (v x ,A(v)v x ) + A(v) (A(v)v x ) x , 
we obtain the modified equation for the Lax-Friedrichs method: 



a, n Ax2 ( Af2 / Air ^ \ 

1H + A(v)v x = — {v xx -^{A (v)v x ) 3 



At 2 ,„ „ \ A 



IX 



^ 2 (DA(v) (A(v)v x ,v x ) - DA(v) (v x , A(v)v x )) j - —I 2 (v). (3.1) 

Note that, in the expression of the modified equations, the only term that depends on 
the choice of the family of paths is the last one: <J> only appears in the expression of ^(w)- 
In a similar way, by a simple (but tedious) calculation we can obtain the equivalent 
equations at third order: 

v t + A(v)v x 

= Ax {jKt Vxx ~ ^Kx Q ^ ~ \ h{ - v ^) + Ax2 (I ( A ( v ) v **)x + \Mv)v xxx 

At 2 

- 2 ((A(v)Q(v)) x DA(v) (Q(v),v x ) - DA(v) (v x ,Q(v))) 



3Ax 2 

At 2 At 2 
J ^D 2 A(v) (A(v)v x , A(v)v x ,v x ) - ^DA(v) (A(v)v x , (A(v)v x ) x ) 

"* {{A{v)I 2 {v)) x + DA{v){h{v),v x )+DA{v){v x ,h{v))) 



w 



4Ax 

--D 2 A(v)(v x ,v x ,v x ) - -h{v) - — 
where we have used the notation 

D 2 A(u)(v 1 ,v 2 ,w) = I J2 vlv 2 m d 2 UkUm A(u) 

\k,m—l 

and 

e(v) = (A 2 (v)v x ) x + DA(v) (A(v)v x , v x ) - DA(v) (v x , A{v)v x ) . 
Furthermore, the expansion of (l2(v)) t is 

{h{v)) t = h,i{v) + h,2{v) + h,3( V )> 
10 



with 



and 



I 2 .i(v)= [ D 2 A(v){A(v)v x ,D ui $-v x ,D ui $ s -v x ) ds 
Jo 

+ / D 2 A(v) (A(v)v x , D Ur <& ■ v x ,D Ur <& s ■ v x ) ds 
Jo 

+ / DA(v) (D 2 UIUI <$> (A(v)v x , v x ) , D ui $ a ■ v x ) ds 
Jo 

+ / DA(v) (D 2 VrV § (A(v)v x ,v x ) , D Ur ^ a ■ v x ) ds, 
Jo 

= [ DA(v) (D Ul <$> ■ (DA(v)(v x , v x ) + A(v)v xx ) , D Ul <S> s ■ v x ) ds 
Jo 

+ / DA(v) (£>„„$ • (DA(v)(v x ,v x ) + A(v)v xx ) , D Ur <S> s ■ v x ) ds 
Jo 

+ / DA(v) (D UI $ • v x , D 2 UiUi <S> s (A(v)v x , v x )) ds 
Jo 

+ [ DA(v) (D Ur $ ■ v x ,D 2 UrU $ s (A(v)v x , v x )) ds, 
Jo 

hM v ) = [ DA ( V ) {D Ul $ ■ v x , D Ul $ s ■ (DA(v){v x , v x ) + A{v)v xx )) ds 
Jo 

+ / DA(v) (£>„ r $ • v x , D Ur $ s ■ (DA(v)(v x , v x ) + A(v)v xx )) ds. 
Jo 

The above formula illustrate the high complexity of the equivalent equation approach 
when dealing with nonconservative schemes. 

3.2. Role of the equivalent equation 

These modified equations are useful to understand why the numerical solutions may 
not converge to the weak solutions of the system. Let us suppose that the family of paths 
<& used in the definition of weak solutions is based on a parabolic regularization of the 
system 

ul + A(u e )ul = e(D(u e )ul) x , (3.2) 
where D is a viscosity matrix which is admissible in the following sense: if ui, u r can be 
connected by a discontinuity satisfying the Rankinc-Hugoniot conditions related to the 
family of paths: 

f 1 d$ 

£(u r -ui)= / A($(s;ui,u r )) — (s;u h u r ) dx (3.3) 



Jo ds 

for some £ € K, then $(s; ui,u r ) is a reparametrization of the solution of the differential 
system: 

- £v' + A(v)v' = (D(v)v')', (3.4) 

with the conditions: 

lim v(£) = U[, lim v(£) = u r . (3.5) 
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Under these hypotheses it can be checked [19] that the function 



u(x,t) = 



ui x < £t, 
u r x > £t. 



(3.6) 



is a weak solution in the sense of Dal Maso, LeFloch, and Murat. 

If now the Lax-Friedrichs scheme is applied to the hyperbolic system, the limits of the 
numerical solutions provide approximations to the vanishing viscosity limits related to 
the regularization (3.1) which is different of (3.2). The difficulty comes from the fact that, 
unlike the conservative case, the vanishing viscosity limits depend on the regularization 
of the problem. Even if, for simplicity, we have only calculated the modified equations 
corresponding to the Lax-Friedrichs scheme, the same difficulty would be present for any 
other scheme involving a numerical viscous term: the numerical solutions approximate 
the vanishing viscosity limit of a modified equation whose regularization terms depend 
both on the chosen family of paths and on the specific form of its viscous terms. 

4. Examples of nonconservative hyperbolic systems 
4.1. A simplified model 



We begin with a hyperbolic system containing nonconservative products, which will be 
used in the following section to perform numerical experiments. We consider the system 



M + Qx = 0, 

„2 s 



qt + 



+ qhh x 



(4.1) 



which has the form w t + A(w)w x — with 

h 



w = 



A(w) = 



1 

i 2 + u 2 h 2u 



and u = q/h. In the region 

SI = {(h, q) | < q, < h< (16g) 1/3 }, 

the system is strictly hyperbolic and all characteristic fields are genuinely nonlinear. The 
eigenvalues of this system are 

Ai = u — h\/u, A2 = u + h\/u, 

and the integral curves of the first and second characteristic fields are 

s/u+ h/2 = const, \Ju — h/2 = const, 

respectively. 

In order to define the jump conditions, the paths connecting the left- and right-hand 
limits w ± = [ft. ± ,(7 ± ] at a shock are chosen to be the union of the segment connecting 
ui~ with w* — [h + ,q~] and the segment connecting w* to w + : 
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<I>(s; w , w + ) = < 



h~ + 2s(h + - h~) 
h+ 

q- + (2s - l)(q+ - q-) 



< s < 1/2; 



1/2 < s < 1; 



(4.2) 



Such paths were introduced in [13] to illustrate certain issues encountered with noncon- 
servative products. The corresponding jump conditions (3.3) are the following: 



m = [q], 



2 



Once the jump conditions have been stated, it is possible to solve the Riemann problem 
for any pair of states which are sufficiently enough. Then, a family of paths satisfying 
(R2) and (R3) is constructed. In Figure 1 the path linking the states wi — [1, 1] T and 
w r — [0.5,0.5] T is depicted. In this case, the solution of the Riemann problem consists 
of a 1-rarefaction connecting wi to an intermediate state w and a 2-shock linking w to 
w r . As a consequence, the path consists of the arc of the 1-integral curve linking wi to w 
and two segments connecting w and w r which are parallel to the axis . The set of states 
that can be connected to wi by an entropy satisfying 1-wave or 2-wave are also depicted. 
In each case Ri denotes i-rarefactions and Si denotes i-shocks. 



4.2. A class of systems of balance laws 



We consider PDE systems of the form: 

w t + F(w) 3 



S(w)a x 



(4.3) 



where the unknown w(x,i) takes values on an open convex set of M. N ; F and 5* are 
regular functions from to R^; and a(x) is a known function from R to K. For an 
account of the existing literature on well-balanced schemes for this class of systems we 
refer to the lecture notes by Bouchut [4], as well as to the recent contribution by Noelle 
et al. [29], and Xin and Shu [33]. 

As pointed out in LeFloch [19] for the Euler equations in nozzle with discontinuous 
cross-section, such a system can be recast in the form of a nonconservative system W t + 
A(W)W X = 0, by introducing 



W 



, A(W) 



J(w) 


-S(w) 









where J(w) denotes the Jacobian matrix of F: 

OF 

J{w) = —(«,). 
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Fig. 1. Path connecting the states w\ and w r (continuous line) and set of states that can be connected 
to wi by an entropy satisfying 1 or 2- wave (dashed lines). 

If J has N different real and non- vanishing eigenvalues Ai(iu), . . . , Xn(w), then the system 
is strictly hyperbolic with eigenvalues 

X 1 (w),...,X N (w),0. 

Clearly the (N + l)-th field is linearly degenerate and, for definiteness, we may assume 
that all other fields are genuinely nonlinear. 

In order to define weak solutions to this nonconservative system, a family of paths 



^(s;W h W r ) 



$ w (s;W h W r ) 
& a {a;W u W r ) 



must be chosen and it is natural to impose the following requirement: if W = [w,cr] T is 
a weak solution to the system and a is a constant, then w must be a weak solution of 
the system of conservation laws: 

w t + F(w) x = 0. 

This requirement is satisfied if the family of paths fulfills the following condition: 
(R4) If Wi and W r are such that 07 = oy = a, then: 



$ a {a;WuW r )=a, a €[0,1]. 



(4.4) 
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In fact any family of paths satisfying (R4) together with the requirements (Rl) and 
(R3) (cf. Section 2.2), leads to the same notion of weak solution for such systems. These 
solutions contain two type of discontinuities : 

- Shock waves across which a is continuous that satisfy the usual Rankine-Hugoniot 
conditions: 

£(to+ - to") = F(w + ) - F(w~). 

- Stationary contact discontinuities placed at the jumps of a and connecting two states 
that belong to the same integral curve of the linearly degenerate field. 

For these systems, the difficulties of convergence commented above are not appreciated 
for the finite difference schemes introduced in Section 2.2. Moreover, the shock waves 
propagating in regions where a is continuous are correctly captured independently of 
the choice of paths. Nevertheless, in order to correctly capture the stationary contact 
discontinuities related to the jumps of a, the numerical schemes have to be based on 
families of paths that satisfy at least the requirement (Rl). Interestingly, this is also the 
requirement necessary to obtain well-balanced schemes. 

The Lax-Friedrichs scheme presented in Section 2.2 cannot be used as it stands for 
systems of balance laws as it does not preserve the equation 

a t = 0. 

In effect, it can be easily verified that the numerical scheme for the variable a reads as 
follows: 

n+1 = a i-l + 

1 2 

In [10] the following modification of the scheme was introduced to get rid of this difficulty: 
instead of (2.9), M™^^ 2 are given by 

M i+l/2 = A i+l/2 ' _ 
1V1 i+l/2 - A i+l/2 \ U i+l 

where 

r ± 1 / . Ax - 



A n ' — - + T n ,4-4" , 

Here, ^4™ +1 / 2 is a Roe linearization and 

Ii+l/2 — ^i+1/2 ' 1^ ' (^r+l/2) ) 

where 3C" +1 y 2 lSi a matrix whose columns are eigenvectors of ^4™ +1 / 2 associated to A™ +1 ^ 2 1 , 

. . . , A" +1 y 2 N , and Id is the diagonal matrix whose j-th coefficient is 1 if A™ +1 ^ 2 ■ ^ 0, or 

if A" +1 / 2 j = 0- Note that if, instead, I™ +1 / 2 1S taken to be equal to the identity matrix, 
then the Lax-Friedrichs scheme presented in Section 2.2 is recovered. 

An important particular example of systems of balance laws is the shallow water system 
governing the flow of a shallow layer of inviscid homogeneous fluid through a straight 
channel with a constant rectangular cross-section: 



dh dq 
dt dx 

dq d ( q 2 q l2 \ , dH 
dt + Yx\T + 2 h =9k dx- 



(4.5) 
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The variable x makes reference to the axis of the channel and t is time; q(x, t) and h(x, t) 
represent the mass-flow and the thickness, respectively; g, the gravity; and H(x), the 
depth measured from a fixed level of reference. 

The eigenvalues of this matrix are Ai = u — c, X 2 = u + c, and 0, where c = \fgh. In 
this case, the equations of the integral curves of the linearly degenerate field are given 
by the equations 

q 2 

q = const, h H - — H = const. (4-6) 

2gh z 

The resonant regime where Ai or A2 vanish are not considered in the present paper; for 
a discussion of the Riemann problem see LeFloch and Thanh [27] . 



4.3. Two-layer shallow water system 



We consider in this paragraph the system of partial differential equations governing 
the one-dimensional flow of two superposed immiscible layers of shallow water fluids over 
a flat bottom topography (see [8] for details) : 

dhi dqi 
dt dx 



0. 



dqi 
dt 
dh 2 
dt 
dqi 
dt 



+ 



+ 



d_ 

dx 

dq 2 

dx 

d_ 

dx 



i. 

hi 



K = -gh 



dh 2 

dx 



(4.7) 



+ 



2 2 



Pi , dhi 

p 2 dx 



In these equations, the index 1 refers to the upper layer and the index 2 to the lower one. 
The fluid is assumed to occupy a straight channel with constant rectangular cross-section 
and constant width. The coordinate x refers to the axis of the channel, t denotes the time 
variable, and g is the gravity. Each layer is assumed to have a constant density pi, i = 1, 2 
(Pi < P2), while the unknowns qt{x,t) and hi(x,t) represent respectively the mass-flow 
and the thickness of the i-th layer at the section of coordinate x at time t. 

System (4.7) can be written in the form (1.1), say w t + A(w) w x — 0, w — w(t,x) £ 
with N = 4 and 



pN 



W 



hi 
h 2 

52 



A(w) 



1 

-u\ + c\ 2m 

rcl 





c\ 
1 

-u\ + c\ 2u 2 



where u, = qijhi represents the averaged velocity of the i-th layer, a — \fg~h~i, i — 1,2, 
and r = ^ . The characteristic equation of the system is 

(A 2 — 2uiA + u\ — ghi) (A 2 — 2M2A + u 2 — gh 2 ) — rg 2 hih 2 . 

Observe that, when r = 0, the eigenvalues are those corresponding to each layer sepa- 
rately. In this situation, the coupling terms do not affect the nature of the system in any 
essential manner. 



1G 



In the case r = 1 (which is the situation arising in many geophysical flows) a first-order 
approximation of the eigenvalues was given in [32] : 



u\h\ + u 2 h 2 



± (g(hi + h 2 )) 



1/2 



(4.8) 



l ext — 



hi + h 2 
u x h 2 + u 2 h x 




h x h 2 



(ui - u 2 ) 2 
g'{hi + h 2 ) 



) 



) 



1/2 




hi + h 2 



{hi + h 2 ) 



(4.9) 



In the former expression, g' = (1 — r)g is the reduced gravity. 

It is not easy to check the genuinely nonlinear character of the 4 characteristic fields, 
as the eigenvalues and eigenvectors can not be written explicitely in a simple manner. 
Nevertheless, this fact is easily proved in the case r = as, in this case, the system 
reduces to two decoupled shallow water systems. As a consequence, using a continuity 
argument, this is also true at least for small values of r. 

From equation (4.9) we can observe that the internal eigenvalues may become complex. 
This situation occurs when they satisfy, approximately, the following inequality: 



In this case, the system loses its hyperbolic character. These situations are related with 
the appearance of shear instabilities that may lead, in real flows, to intense mixing of the 
two layers. While, in practice, this mixture partially dissipates the energy, in numerical 
experiments these interface disturbances grow and overwhelm the solution. Clearly, we 
cannot expect to simulate these phenomena with a two-immisciblc-layer model. There- 
fore, the above inequality in fact gives the range of validity of a model based on the 
equations (4.7), if viscosity effects are neglected. In this work only the case where the 
matrix A(w) has real eigenvalues is considered, i.e. the system is supposed to be strictly 
hyperbolic. 

The jump conditions (3.3) related to the choice of a family of paths 



(mi - u 2 ) 2 
cfihr + h 2 ) 



> 1. 



$(s; wi, w r ) = 



$ hl (s;wi,w r ) 
® qi (s;wi,w r ) 

$ q2 (s;w h w r ) 



read 
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m-h l 1 )=q{-q[, 

+9 ( ^h 1 {s;w h w r )^ }2L {s;wuWr) ds, 
Jo ds 



ail ~ <Q 



is {<&) (fir , 5/f,rs2 S/.Jn2 



/'2 



+gr I $ h2 (s;wi,w r )—— J -(s;w l ,w r )ds. 



(4.10) 



Observe that these conditions are independent of the choice of & qi (s : wi,w r ), i = 1,2. 




1.2 13 1.4 1.5 1.6 1.7 18 1.5 



Fig. 2. Test case 5.1. Hugoniot curves: exact (continuous red line) and numerical (line with dots) obtained 
with Roc scheme 



5. Numerical experiments 

5.1. A simplified system 

We will demonstrate here that, in general, a difference scheme for a general noncon- 
servative hyperbolic system (1.1) does not converge to the exact solution u, that is, with 
the notation already introduced in previous sections we claim that 

v = lim u A ^ u. (5.1) 
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This fact was first observed for nonconservative schemes for scalar equations in Hou 
and LeFloch [16] and, in the context of nonclassical shocks generated by diffusion and 
dispersion, in Hayes and LeFloch [15,25]. Here, we observe (5.1) for finite difference 
approximations of nonconservative systems. 



— Roe 
Gosynov 




1.2 1.3 1.4 1.5 



Fig. 3. Test case 5.1. Hugoniot curves: exact (continuous red line), Roe (blue line with squares) and 
Godunov (green line with circles) 



We begin with the convergence of Roe and Godunov methods for (4.1). The Roe 
method considered here is consistent with the family of paths given by (4.2) for every 
pair of states. The corresponding Roe matrix is as follows: 



4™ — 

A i+l 2 - 







~( U "+l/2) 2 + l"K+l/2 2u "+l/2 



where 



1 



H+l/2 



i+i/2 — ^(K + hi+l)- 



On the other hand, the Godunov method considered here is based on the family of 
paths described in Section 4.1: it satisfies (R2)-(R3) and coincides with (4.2) for pair of 
states that can be linked by a shock. 

We consider the Rankine-Hugoniot curve composed by the states w r that can connected 
with wi = [1; 1] T by a 1-shock, which is given by: 



1 



1 



2h r 



(K - 1) 



(5.2) 



Figure 2 shows a plot of the curve (5.2) in the plane h — q (continuous red line). 
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(b) h (zoom) at t = 0.5. 



(a) h at t = 0.5 



Fig. 4. Test case 5.1. Solution of Riemann problem (h) (5.3) at time t = 0.5: Exact (continuous red line), 
Godunov (green line with circles) and Glimm (blue line with squares) 




(a) q at t = 0.5 (b) q (zoom) at t = 0.5. 

Fig. 5. Test case 5.1. Solution of Riemann problem (q) (5.3) at time t = 0.5: Exact (continuous red line), 
Godunov (green line with circles) and Glimm (blue line with squares) 

Next, Roe method is used to solve numerically a family of Riemann problems in which 
the left state is wi and the right state w r runs on the Hugoniot curve (5.2). The speed of 
propagation and the limit states of the 1-shock are computed in the numerical solution 
by using the first divided difference as a smooth indicator. This computation has been 
performed using four meshes with decreasing mesh step (Ax — 0.002, Ax — 0.001, Ax = 
0.0005 and Ax — 0.00025). The CFL parameter is set to 0.9. The numerical Hugoniot 
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curves obtained in this way (line with dots) are compared with the exact one (continuous 
red line) in Figure 2. It can be observed that the numerical Hugoniot curves converge, 
but the limit is not the exact one. 

The same behavior is observed for Godunov method: the numerical Hugoniot-curves 
converge but the limit is not the exact one: in Figure 3 the exact Hugoniot curve (con- 
tinuous red line) is compared with those computed with Roe (blue line with squares) 
and Godunov (green line with circles) methods with Ax = 0.001. The CLF parameter is 
set to 0.9 for Roe and 0.5 for Godunov. This choice ensures that the Godunov method 
corresponds to advance in time by exactly solving the Riemann problems and taking the 
averages of the solutions at the cells. 

Finally, we compare Godunov and Glimm methods. We consider a Riemann problem 
with initial conditions 




x < 0; 

(5.3) 




w(x, 0) 



where q r is given by (5.2) for h r = 1.8 (q r = 0.530039370688997). The exact solution 
consists thus of a 1-shock linking the states. In Figures 4 and 5 the exact solution at time 
t = 0.5 is compared with numerical solutions obtained with Godunov (green line with 
circles) and Glimm (blue line with squares) methods with Ax = 0.001 and CFL=0.5. 
Note how Godunov method introduces a 2-rarefaction in the computed solution. 

Observe that the first equation of (4.1) is a conservation law. According to this equa- 
tion, for A > large enough, the exact solution of the Riemann problem has to satisfy 
the following conservation property: 

rA r A 



h(x,t)dx= / h(x, 0) dx + t(l — q r ). 

-A J -A 



In Figure 6 we investigate the conservation property of both Godunov and Glimm meth- 
ods: we fix A and compare the exact value of the integral of h at time t n with its numerical 
approximations. Note that Godunov method satisfies the conservation property exactly. 



5.2. Shallow water system 

In this section we consider the discretization of the shallow water system by means 
of Roe and the modified Lax-Fricdrichs schemes applied to the formulation (1.1) of the 
problem. In the first test, we check that, for continuous bottom functions H, the shock 
waves are correctly captured for the schemes even if the numerical schemes are formally 
consistent with a simple family of paths. In the second one, stationary contact discontinu- 
ities placed at the jumps of H are considered: we check that they are correctly captured 
only when the family of paths satisfies the property (Rl). 
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0.1 0.2 0.3 0.4 0.5 



Fig. 6. Test case 5.1. Evolution of the 'mass': exact (continuous red line), Glimm (blue line with squares) 
and Godunov (green line with circles) 

5.2.1. Dam-break problem over a non-flat bottom topography 

The axis of the channel is the interval [0, 10] and the bottom topography is given by 
the function H(x) = 1 — 0.5e~ < - ;E ~ 5 - ) . The initial condition is q — and 

{H(x), x>4; 
H(x)+0.5, x<4. 

The final time is t = 0.6. Free boundary conditions are considered. The CFL parameter 
is set to 0.9. 

We consider a Roe scheme and a modified Lax-Friedrichs scheme which are consistent 
with the family of segments. Figure 7 shows the bottom topography and the free surface 
computed for both schemes using three meshes with increasing number of cells (800, 1600 
and 3200 cells respectively). Both schemes converge to the same solution. Moreover, the 
speed of propagation £ and the limit states w~ and w + of the shock in the numerical 
solutions have been computed for both schemes by using a fine mesh of 32000 cells and 
the first divided difference as a smooth indicator. The value of the residual \£,(w + — w~) — 
F(w + ) + F(w~)\ obtained for the well-balanced Lax-Friedrichs scheme is 0.008 and for 
Roe scheme, 0.006. 



5.2.2. Stationary contact discontinuities 

In this test we study the approximation of stationary contact discontinuities. Following 
the discussion in Section 4.2, such a discontinuity has to connect two states belonging 
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(a) Bottom topography and free surface at t=0.6. 



(b) Free surface (zoom) at t=0.6. 




Fig. 7. Test case 5.2.1. Dam-break problem: bottom topography and free surface at t=0.6. 

to a same curve of the family (4.6). As the family of segments only satisfies (Rl) for the 
particular case of the curves 

q = 0, h — H = const, 

the Roe and the modified Lax-Friedrichs schemes based on the family of segments arc only 
expected to capture correctly stationary contact discontinuities corresponding to water 
at rest over a discontinuous bottom. To check this in practice, we consider a channel 
whose axis is the interval [—5,5] and whose bottom is given by the function 



H(x) 

We consider the initial condition: 

Wi 



where h r has been calculated so that both states belong to the same integral curve 
(4.6) {h r = 0.7892441190408083). The exact solution of this Riemann problem is thus a 
stationary contact discontinuity. 

We have applied both schemes to this Riemann problem. As boundary condition, 
the state wi is imposed upstream and free boundary conditions downstream. The CLF 
parameter is set to 0.9. Figures 8 and 9 show the stationary solutions obtained with 
both schemes using three meshes with increasing number of cells (100, 200 and 400 cells 
respectively). As expected, the numerical solutions do not converge to the exact solution. 
Note that both schemes converge to the same discontinuous function. 



h 




l 


q 




V*9 



Let us check what happens if the numerical schemes are based on a family of paths 
satisfying (Rl) for every integral curve of the linearly degenerate field. We consider the 
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(a) Bottom topography and free surface (stationary (b) Bottom topography and free surface (station- 
solution), ary solution): zoom 

Fig. 8. Test case 5.2.2. Stationary contact discontinuity: Comparison between the modified Lax-Friedrichs 
and Roe schemes and the exact solution (free surface). 





Fig. 9. Test case 5.2.2. Stationary contact discontinuity: Comparison between the well-balanced 
Lax-Friedrichs and Roe schemes and the exact solution (discharge). 

family of paths constructed in [9] in order to design a Roe scheme which is well-balanced 
for every smooth stationary solution. The path connecting two states Wi and W r consists 
of an arc of one of the curves of the family (4.6) and a segment lying on a plane H — 
const: see [9] for details. Once the corresponding Roe matrix has been calculated, the 
construction of the modified Lax-Friedrichs scheme is straightforward. 

Figure 10 shows the comparison between Roe, the modified Lax-Friedrichs scheme, and 
the exact solution for three meshes with increasing number of cells (100, 200 and 400 
cells respectively). As expected, the stationary contact discontinuity is exactly captured. 
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(a) Bottom topography and free surface (stationary (b) Bottom topography and free surface (station- 
solution), ary solution): zoom 

Fig. 10. Test case 5.2.2. Stationary contact discontinuity: Comparison between the exactly well-balanced 
Lax-Friedrichs and Roe schemes and the exact solution (free surface). 

5.3. Two-layer shallow water system 

5.3.1. Approximation of internal shocks 

In this section we consider the discretization of the homogeneous two-layer shallow 
water system (that is, H = est.) by means of Roe and Lax-Friedrichs schemes. We check 
that the numerical solutions do not converge to the weak solutions involving shocks even 
when the same family of paths is used both for the definition of the jump conditions and 
the construction of the numerical scheme. 

We begin with the family of segments for the definition of the jump condition (4.10). 
In this case: 

9 hl (s;wi,w r )—Q—(s;w[,w r ) as = (h 2 - h 2 ), 

, . d$ hl , s , hi + h\ ., _ , i. 

9 h2 (s;wi,w r )— — (s;wi,w r ) as = - ftj. 

os 2 

In all of the cases considered here the order of the eigenvalues of the system is: 
Moreover 

« \Xtt\> (5-4) 
We consider first a Lax-Friedrichs scheme consistent with the family of segments. Some 

easy calculations show that, for the particular choice of the family of segments, the last 

term in the modified equation (3.1) vanishes. 

The goal is to compare the exact and the numerical Hugoniot curves corresponding to 

one of the internal characteristic fields; i.e. the fields related to the eigenvalues Xf nt ). We 

proceed as follows: the state 
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K 




0.392034161025472 






-0.198826959396196 


hi 




1.588829011097482 






0.186046955388750 



is fixed. Then, we compute the Hugoniot curve corresponding to the 'left 'states wi that 
can be connected with w r with a 3-shock. To do this, we use the speed of the shock £ as 
a parameter and, for each value of £ we solve the non-linear system (4.10). In Figure 11 
we show the projection (continuous blue line) of the computed Hugoniot curve onto the 
planes h\ — q\ (left) and hi — q-i (right), respectively. 

Next, we solve numerically a family of Riemann problems in which the right state is 
w r , while wi runs on the Hugoniot curve. Using the first divided difference as a smooth 
indicator, the speed of propagation and the limit states of the shock corresponding to 
the eigenvalue Xf nt is determined in the numerical solution. These calculations have been 
performed by using four meshes with decreasing step ( Ax = 0.002, 0.001, 0.0005 and 
0.00025). The numerical Hugoniot curves so obtained are compared with the exact one 
in Figure 11. Observe that the numerical solutions converge, but the limit is not a weak 
solution according to the chosen family of paths. Nevertheless, if w r and wi are close 
enough, both curves are very close. The same behaviour can be observed if the shock 
speed is close to zero: this situation corresponds in the figure to the intersections of the 
curves. 




(a) Hugoniot curves (projection onto the plane h\ — (b) Hugoniot curves (projection onto the plane 
qi): exact (continuous blue line) and numerical h^ — q^)'- exact (continuous blue line) and numerical 
(lines with dots) (lines with dots) 

Fig. 11. Test case 5.3.1. Hugoniot curves: exact (continuous blue line) and numerical (lines with dots). 

A similar behavior is observed for Roe scheme: the numerical approximations converge 
but the limit is not a weak solution according to the family of segments. Nevertheless, 
as the numerical viscosity of the scheme is smaller than that corresponding to Lax- 
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Friedrichs, the results are expected to be closer to the exact solution: in Figure 12 we 
compare the exact Hugoniot curve with those computed with Lax-Friedrichs and Roe 
schemes using a mesh with 10000 cells. 





(a) Hugoniot curves (projection onto the plane hi — (b) Hugoniot curves (projection onto the plane hi — 
qi): exact (continuous red line) and numerical (lines qi): exact (continuous red line) and numerical (lines 
with dots) with dots) 



Fig. 12. Test case 5.3.1. Hugoniot curves: exact (continuous red line) and numerical (lines with dots). 



5.3.2. Approximation of external shocks 

Due to the inequality (5.4)the CFL condition adjusts the numerical velocity to the 
external eigenvalues. As a consequence the effects of the numerical viscosity are much 
stronger for internal shocks and thus external shocks are expected to be better captured 
with Lax-Friedrichs or Roe schemes. 

In order to check this, we proceed as in the previous test case: the state w r given by 



0.257381469591567 
0.444901654188681 
0.110306344093418 
0.190672137450279 



w r 



is fixed, and then the Hugoniot curves corresponding to the 'left 'states (wi) that can be 
connected with w r with a shock related to the eigenvalue X~ xt are computed by solving 
the non- linear system (4.10) using £ as a parameter. 

In Figure 13 we compare the exact Hugoniot curve with those computed with Lax- 
Friedrichs and Roe scheme using a mesh with 2000 cells. Note how the curves are now 
much closer to each other than they were in the previous test case. 
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(a) Hugoniot curves (projection onto the plane hi — (b) Hugoniot curves (zoom of the projection onto 
qi): exact (continuous red line) and numerical (lines the plane h± — qi): exact (continuous red line) and 
with dots) numerical (lines with dots) 




(c) Hugoniot curves (projection onto the plane hi — (d) Hugoniot curves (zoom of the projection onto 
52): exact (continuous red line) and numerical (lines the plane /12 — 92): exact (continuous red line) and 
with dots) numerical (lines with dots) 

Fig. 13. Test case 5.3.2. Hugoniot curves for Lax-Friedrichs and Roe schemes: exact (continuous red line) 
and numerical (lines with dots). 

5.4. Influence of the family of paths 

In this test we study the influence of small changes in the family of paths both in the 
weak and the numerical solutions. We consider now the jump conditions related to a 
different family of paths Q e h (s; wi, w r ), i = 1,2. These curves are chosen so that: 

hi = ®h 1 ( s '> w h w r), se[0,l]; 
hi = $l 2 (s;w u w r ), se[0,l]; 



28 



is a parameterization of the curve 



XKf-(h\f 



h r 2 -h l 2 



h\-h[ {hiY-{h{f 

The jump conditions are now (4.10), with 

'U dq>h A( W (^j) 2 (3 + 46) + 2(3 + 2e)h[h\ + (3 + Ae)(h\f 

d ® h A( w h\((3 + 4e)h 2 + (3 + 2e)h r 2 )+h[((3 + 2e)h 2 + (3^ 



Ae)h\) 



Observe that the jump conditions of the previous tests are recovered for e = 0. 

As in Section 5.3.1, the state w r given by (5.5) is fixed and the exact Hugoniot curves 
are computed for different values of e. In Figure 14 the projections of the exact Hugoniot 
curves onto the planes hi — q± (left) and h 2 — 92 (right) for the values e = 0.00, 0.01, 0.02, 
0.03, 0.04, 0.05, are shown. 





(a) Hugoniot curves (projection onto the plane hi - 



(b) Hugoniot curves (projection onto the plane hs- 
12) 



Fig. 14. Test case 5.4. Hugoniot curves: e e {0.0, 0.01, 0.02, 0.03, 0.04, 0.05}. 



Next, we consider the Lax-Friedrichs schemes related to the new choice of the family 
of paths. An easy calculation shows again that, for this new family of paths, the last 
term of the modified equation (3.1) also vanishes. As a consequence, the second order 
modified equation is independent of e: it coincides with the one corresponding to the 
family of segments. 

We proceed as in the previous test case to compute the Hugoniot curves corresponding 
to the numerical scheme. In Figure 15 the curves obtained using a mesh with step Arc = 
0.001 for e = 0.00, 0.01, 0.02, 0.03, 0.04, 0.05, are depicted. Observe that all of the curves 
coincide: the Hugoniot curves obtained with the Lax-Friedrichs scheme for the different 
values of e are reparameterizations of the same curve. This fact agrees with the fact that 
the second order modified equation is independent of e. 
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Fig. 15. Test case 5.4. Hugoniot curves for Lax-Friedrichs scheme with 
e 6 {0.0, 0.01, 0.02, 0.03, 0.04, 0.05}: exact (red lines) and numerical (lines with dots). 



Figure 16 shows the projections of the Hugoniot curve corresponding to Roe scheme 
for e = 0.00, 0.01, 0.02, 0.03, 0.04, 0.05 (line with dots) onto the planes hi - qi and 
hi — q2- Note that, in this case, the obtained curves depend on e, but again the numerical 
solutions do not converge to the weak solutions corresponding to the chosen family of 
paths. Nevertheless, as in the previous test case, if w r and Wi are close enough or the 
shock speed is close to zero, the exact and the numerical Hugoniot curves are also close. 
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(a) 3-shock hugoniot curves (projection onto the (b) 3-shok hugoniot curves (projection onto the 
plane hi — qi) plane hi — 92) 



Fig. 16. Test case 5.4. 3-shock hugoniot curves (zoom) for Roe scheme with 
e 6 {0.0, 0.01, 0.02, 0.03, 0.04, 0.05}: exact (red lines) and numerical (lines with dots). 

6. Conclusions 

When a hyperbolic system with nonconservative products and genuinely nonlinear 
fields is discretized, in order to be sure that the numerical approximations converge to 
a function which is a classical solution where it is smooth and whose discontinuities are 
in good agreement with the physics of the problem, the following steps should be taken 
[19,16]: 

- First, choose a regularization of the system which is consistent with the physics of the 
problem. 

- Next, determine the DLM family of paths consistent with this regularization. 

- Finally, design a numerical scheme whose solutions converge to weak solutions associ- 
ated with this family of paths. 

In practice, this strategy may be difficult to follow, since the actual calculation of a 
family of paths requires calculating regularized shock profiles associated with the given 
regularization. On the other hand, the convergence of the numerical solutions to the 
correct weak solutions is known for the Glimm scheme and the front tracking method 
[24], only; the implementation of these methods can be time consuming since they require 
the explicit knowledge of the corresponding Riemann solver. 

In fact, when the nonconservative model under consideration is a simplified version of 
a more complex (but conservative) model — as is the case for the two-layer shallow water 
system — the above strategy may end up being more costly than than solving directly 
the more complex one. In these cases, the use of a numerical strategy based on a direct 
discretization of the nonconservative system by means of finite difference scheme which 
is formally path-consistent is advisable and may have the following advantages: 

- The numerical solutions is formally consistent with the definition of the nonconser- 
vative product in the sense of Dal Maso, LeFloch, and Murat and, in turn, in the 
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special case that the system admits a conservative subsystem, the numerical scheme 
is conservative for that subsystem in the sense of Lax. 

- The approximations of the shocks provided by the schemes are consistent with a reg- 
ularization of the system with higher-order terms that vanish as Ax tends to 0. (Ob- 
viously, the main drawback is that this rcgularization depends on the chosen family 
of paths and on the numerical scheme itself. This is issue dealt with in the present 
paper.) 

- As originally pointed out by Hou and LeFloch [16] in the (simpler) case of scalar 
hyperbolic equations, the convergence error, measured in terms of our convergence 
error measure or in terms of the Hugoniot curves, is noticeable for very fine meshes, 
for discontinuities of great amplitude, and/or for large-time simulations, only. 

This strategy is extendable to high-order methods or to multidimensional problems, 
as developed, together with collaborators, by Coquel [1,2,11] and Pares [5,6]. 
The convergence error should also be compared with the experimental error. In the 
case of the two-layer shallow water system, the shocks captured by Roe scheme and the 
family of straightlincs have been found [7] to be in good agreement with the experimental 
measurements of internal bores in the Strait of Gibraltar, despite of the simplicity of the 
family of paths. 

In certain special situations, the convergence error measure is found to vanish iden- 
tically. This is the case of systems whose nonconservative product is associated with a 
linearly degenerate field: for schemes that arc formally consistent with a family of paths 
satisfying the condition (Rl), then all of the discontinuities are correctly approximated 
and the scheme does converge to exact solutions. The discussion of linearly degenerate 
fields associated with nonconservative products was discussed earlier in [19,26] from the 
theoretical standpoint and in [4] from the numerical standpoint. This problem may also 
exhibit an additional difficulty, the resonance problem, if one of the eigenvalues of the 
Jacobian matrix vanishes, and weak solutions may not be uniquely determined by their 
initial data, so that the limiting numerical solutions may depend both on the family of 
paths and the numerical scheme itself. 
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